function [P_S, F_es, wprimeS, educ_mat] = choiceProb(Ss, rs)
    global  a bet0 bet1 pe tau w ub lb N_w N_a sigmae b ...
            thetas w_bar unemp xi eta  gamma share curve w_gamma idx_pe
    
    % Pre-allocate matrices
    N_n      = length(Ss);
    u_s      = zeros(N_w, N_a, length(Ss), 2);
    educ_mat = zeros(N_w, N_a, length(Ss));
    F_es     = zeros(N_w, N_a, length(Ss));
    wprimeS  = zeros(N_w, N_a, length(Ss));
    [~,id]   = min(rs);
    pe       = zeros(N_w,1);
  
    % Solve education decision for each potential destination neighborhood 
    for c = 1:length(Ss)
        % Return to unit of education
        E_S = eta * (a' .* (bet0 + bet1*Ss(c)).^xi)  .* f_w(w) ; 
        % Wage additive constant
        btilde = unemp + b .* f_w(w) + pe .* E_S ;
        
        % Solve for optimal education
        % Quadratic Solution as initial guess
        discrim = btilde.^2 + (3./tau)*((E_S).^2).*(w - rs(c));
        e_init = max(min((-btilde + sqrt(max(discrim,0)))./(3 .* (E_S)) .* (discrim > 0), ub), lb) + 0.1;
      
       % Newtons Method
        for iter = 1:230
           % Precompute
           evar = e_init.^(gamma - 1);
           evar2 = e_init.^gamma;
           %Residual
           foc =  gamma*tau* evar .* btilde + gamma*E_S.*tau .* evar2 - (w - rs(c) - tau*evar2) .* E_S * w_gamma; 
           foc = foc .* (w - rs(c) > 0);
           foc(isnan(foc)) = 0;
           % Residual derivative
           dfoc =  gamma*(gamma - 1)*tau*e_init.^(gamma - 2) .* btilde + gamma^2*E_S.*tau .* evar  + gamma*tau*evar.* E_S * w_gamma; 
           e_init = (e_init - foc ./ dfoc) .* (w - rs(c) > 0);
           e_init(isnan(e_init )) = 0;
           % Check solution
           if max(abs(foc(:))) < 5e-4
               break
           elseif iter == 130
               display(max(abs(foc(:))))
               error("Did not converge in education")
           end
        end
   
        e_S = e_init;
        e_S = min(max(e_S, lb), ub);
        
        % For use in other programs
        F_es(:,:,c)     = (pe + e_S) .* E_S ./  f_w(w) ;
        educ_mat(:,:,c) = pe + e_S;
        wprimeS(:,:,c)  = (unemp + (b  + F_es(:,:,c)) .* f_w(w)) ;
        u_s(:,:, c, 1)  = (max((w - rs(c) - tau*e_S.^gamma).*( wprimeS(:,:,c).^w_gamma ), 0) );
        
        % If wage is too low to afford housing, allocate to cheapest neighborhood
        if (c == length(Ss))
            mysub = sum(u_s, 3:4);
            [mrow, mcol] = ind2sub(size(mysub), find(mysub == 0));
            u_s(mrow,mcol, id,1) = 1e-20;
        end
    end

    % Amenities shock
    for c = 1:length(Ss)
        u_s(:,:,c,2) = u_s(:,:,c,1) .* exp(log(thetas(c)));
    end

    % Calculate choice probabilities
    my_log_diff = (log(u_s) - max(log(u_s), [], 3)) ; 
    numer = exp(my_log_diff .* sigmae ); 
    P_S = numer ./sum(numer , 3);
      
end